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Abstract: 



A variational method for computing conformational properties of molecules with Lennard-Jones potentials 
for the monomer-monomer interactions is presented. The approach is tailored to deal with angular degrees 
of freedom, rotors, and consists in the iterative solution of a set of deterministic equations with annealing 
in temperature. 

The singular short-distance behaviour of the Lennard-Jones potential is adiabatically switched on in order 
to obtain stable convergence. 

As testbeds for the approach two distinct ensembles of molecules are used, characterized by a roughly 
dense-packed ore a more elongated ground state. For the latter, problems are generated from natural 
frequencies of occurrence of amino acids and phenomenologically determined potential parameters; they 
seem to represent less disorder than was previously assumed in synthetic protein studies. 

For the dense-packed problems in particular, the variational algorithm clearly outperforms a gradient 
descent method in terms of minimal energies. Although it cannot compete with a careful simulating 
annealing algorithm, the variational approach requires only a tiny fraction of the computer time. 

Issues and results when applying the method to polyelectrolytes at a finite temperature are also briefly 
discussed. 
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1 Introduction 



The determination of configurations of long molecular chains often is a difficult problem. 
For polyelectrolytes, consisting of identical charged monomers interacting with Coulomb 
repulsion forces, the ground state is trivial - the monomers form a straight line. The 
challenge here lies in predicting statistical quantities for finite temperature configurations; 
thus, the thermodynamics of the system is crucial. For proteins, modelled as a sequence of 
point-like amino acids interacting with effective pair potentials having local minima, the 
situation is somewhat different. Here the groundstate is nontrivial; the energy landscape 
is typically plagued with many local minima. This is the main difficulty here, while the 
finite temperature properties often are just considered as minor perturbations around the 
ground state. 

Optimization problems with many local minima as in the the protein case are notoriously 
difficult, and elaborate methods to search the phase space efficiently have been devised. 
One such method is Simulated Annealing (SA) [1], where noise is introduced to emulate 
a Boltzmann distribution; this enables the system to escape from local minima. Unfor- 
tunately, this procedure is quite CPU demanding. The formal temperature in such an 
approach is merely an artificial parameter; it need not be identified with the physical 
temperature of the system. 

For some optimization problems it has turned out profitable to abandon stochastic methods 
like SA in favour of deterministic approaches based on the iterative solution of equations 
originating from a variational scheme. Successful results with such approaches were re- 
ported for combinatorial optimization problems that can be mapped onto spin systems 
[2, 3]. 

Variational methods have recently also been used to compute configurational properties 
of polyelectrolytes. In refs. [4, 5], harmonic trial potentials were used with widths and 
positions as variational parameters. This approach works well as long as the potentials to 
be approximated have a milder divergence than 1/r 3 . In the Coulomb case (screened or 
unscreened), which was treated in refs. [4, 5], this requirement is fulfilled. However, for 
the strongly diverging (1/r 12 ) Lennard-Jones (LJ) potential, occuring in protein models, 
a harmonic Ansatz will lead to divergent integrals, as will indeed any smooth Ansatz. 

This problem could be overcome in two different ways: In principle, one could use a modified 
Ansatz for the trial potential, that leads to finite integrals. However, this is difficult to 
achieve while maintaining the computational simplicity needed for a competitive algorithm. 
Alternatively, one could keep the harmonic Ansatz, and instead modify the LJ-potential 
to make it less singular. 

The approach of refs [4, 5] applies only to molecules with flexible bonds, however; in this 
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Figure 1: (a). Elementary moves on the unit sphere, (b). Evolution of a MF rotor initialized close to the 
center, for T < T c . 



paper we will consider a slightly simplified model with rigid bonds, and the dynamics 
hence limited to the angular degrees of freedom. This simplification is motivated by the 
assumption that the flexibility of bonds is of minor importance for obtaining an accurate 
spatial structure. Indeed, several calculations using real chains of amino acids focus entirely 
on the angular degrees of freedom [6] . 

Thus, we are facing a minimization problem in terms of an energy-function of angular 
variables. These will be represented by a set of unit vectors, rotors, which can be seen 
as continuous generalizations of discrete Ising spins. In algorithms like gradient descent 
(GD) and SA elementary moves are made with the rotors on-shell, ie. restricted to the 
unit sphere (see fig. la). 

In ref. [10] the variational mean-field approach, commonly used for discrete spin-systems, 
was generalized by introducing mean field rotors; these can explore the interior of the sphere 
and can be interpreted as thermal averages of on-shell rotors. The method was successfully 
explored on the minimal energy problem for charges on a sphere. The corresponding mean 
field equations are iteratively solved as the temperature T is lowered. Above a critical 
temperature T c the system relaxes to a fixed point with all rotors at the center of the 
sphere. As the temperature is lowered the rotors approach the "shell" (see fig. lb). Thus, 
in the variational approach, rather than fully or partly exploring the configuration space, 
the variables "feel" their ways off-shell in a fuzzy manner towards good on-shell solutions. 

The main purpose of this paper is to apply a similar technique to proteins, modelled as a 
chain of monomers connected by rigid bonds and interacting via LJ-potentials, with the 
aim of finding the ground state. 

Any updating algorithm faces problems with the steepness of the LJ-potential at short 
distances. We have developed an adiabatic regularization procedure that efficiently handles 
the short-distance problem, which is useful for any method, not just the variational rotor 
approach. 
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Computational simplicity is gained, at the price of sacrificing uniqueness and physical 
interpretability at non-zero T, by approximating expectation values according to 



(1) 



The resulting algorithm is shown to yield a GD algorithm in the T — > limit. 

As LJ-potential testbeds we do not consider real- world proteins. For the purpose of algo- 
rithmic development and studies it suffices to study synthetic systems. We have chosen to 
study two extremes, by considering systems with a dense-packed or elongated ground state, 
respectively. The latter were generated from empirical independent amino-acid probabili- 
ties and phenomenologically determined pairwise forces. As a by-product, it turns out that 
the resulting LJ-parameters give rise to far less disordered systems than has been assumed 
in some recent generic investigations of spectra and stability issues [8, 9]. 

For these testbeds the variational rotor approach is explored for problem sizes N ranging 
from 10 to 40. For the dense-packed model the variational approach compares favourably 
with GD with respect to reaching low energy states, while for the elongated model the 
corresponding gain is smaller. For both models the CPU time consumption is much lower 
than for SA. 

For completeness, we apply the variational rotor approach also to a polyelectrolyte (Coulomb) 
problem at T ^ 0, and discuss the limitations of the approach here. 

This paper is organized as follows. In Sect. 2 the variational (mean field) formalism is 
derived. The generation of synthetic proteins and their LJ-potential couplings are discussed 
in Sect. 3. The regularization of the LJ-potential is treated in Sect. 4. In Sect. 5 we 
present the numerical procedures and results for the LJ-potential, and Sect. 6 contains the 
polyelectrolyte application. A brief summary and outlook can be found in Sect. 7. 



2 The Variational Approach 
2.1 Proper Variational Approach 

Limiting ourselves to angular degrees of freedom, we consider an energy function 



to be minimized w.r.t. the directions of a set of N distinct D- dimensional unit vectors S; 



E = E(s 1 ,...,s N ) 



(2) 



(rotors) 



(3) 
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In instances where the energy landscape contains many local minima, one would typically 
employ a stochastic technique like SA. In [10] a mean field method was developed and 
numerically explored for the problem of placing charges on a sphere. We will here generalize 
this technique, and apply it to the case of a protein model with LJ pair-potentials. 

For spin systems the mean field approximation can be derived in (at least) three concep- 
tionally distinct ways; from a variational principle, from a saddlepoint approximation, or 
using an intuitive physical argument. We will here briefly discuss the variational derivation; 
for the saddle-point approach we refer the reader to e.g. ref. [10]. 

The variational approach is based on an effective energy Ansatz, Ey , which is linear in the 
spins 

E v (s 1 , . . . ,s N ) = -J2 u i ■ s i (4) 

i 

where the (real, unconstrained) coefficient vectors U; are to be considered as variational 
parameters. 

Based on the corresponding Boltzmann distribution (oc exp(— Ey/T)), a free energy, Fy, 
can be defined, 

Fy(u 1 ,...,U N ) = (E)y-ST, (5) 

where S is the entropy. Fy is bounded from below by the true equilibrium free energy, F, 
based on the proper Boltzmann distribution oc exp(— E/T). 

The variational free energy Fy can be written as 

Fy = -Tlog Zy + (E — Ey)y (6) 

where Zy is the variational partition function, while ()y refers to averages w.r.t. the 
variational Boltzmann distribution. By minimizing Fy with respect to the parameters U;, 
the variational equations result. These are best expressed in terms of the mean fields 

V; = (si) v (7) 

which approximate the exact averages (s;). The mean fields are simple functions of the 
coefficients U;, 

v i = ^(|u i |/r) = g (u i /r) (8) 

i ^* i 

where for the case of three dimensions g is given by g(x) = coth(a;) — 1/x; its graph is shown 
in fig. 2. Note that g(0) = 0, and that g(x) — > 1, when x — > oo; this generic qualitative 
behaviour holds for any number of dimensions, and implies in particular that as T — > the 
MF rotors {v;} go on-shell, and can be identified with a low-energy configuration {s;}. 

The variation of the entropy term ST yields 

d(ST) = -J2"idvi (9) 



4 



1 



O 2 4 6 8 10 

Figure 2: The function g(x) for the case of three dimensions. 



and we obtain the variational equations in the form 

Vi = &{-V Vi (E) v /T) 



(10) 



2.2 Modified Variational Approach 



For a strongly singular potential, like LJ, things are complicated by the fact that {E)y 
is a diverging integral. In other cases, the corresponding integral may be convergent, but 
difficult to evaluate. These difficulties can be remedied by making the apparently crude 
replacement 

{E(s 1 ,...,s N )) v -> E({s 1 ) v ,...,{s N )v) = E(y u ...,v N ) (11) 

in the expression for Fy. This is justified as long as we are only interested in the ground- 
state, which dominates for T — > where the fluctuations vanish and the above approxima- 
tion becomes exact. For a finite T, the replacement in eq. (11) is more questionable, and 
its use must be justified by other means. 

Minimization of the thus modified Fy with respect to the trial parameters U; (or Vj) yields 
a modified set of equations: 

v i = g(-V i E(v 1 ,...,v N )/T) (12) 
In principle, these could be iterated to find a (local) minimum of Fy. 

However, a question of uniqueness arises: since E really is defined only on-shell, i.e. for 
|sj| = 1, off-shell values are not uniquely defined. By adding eg. a term oc v? — 1 which 
is zero on-shell, we can alter the off-shell behaviour of E, and hence the solutions to eq. 
(12). Thus, these make sense only in the T — > limit, where the V; are forced on-shell, 
due to the asymptotic behaviour of g(). 

Another problem is a possible lack of stability of the iterative dynamics of eq. (12). This 
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can be remedied by formally adding a stabilizing term to the energy of the form —j3 /2 J2i "vf ) 
which is just a constant on-shell. 

We then obtain a modified set of variational equations: 



to be iterated. For T — > this turns into a kind of an (on-shell) gradient descent with /3 
corresponding to a reciprocal step size. This of course could have been obtained in a much 
simpler way. However, we aim at an annealing approach: start iterating with a high T, 
for which the modified Fy typically is minimized by a fixpoint with V; 0. Then keep 
iterating while slowly letting T — > 0. The value of (3 should be chosen to stabilize the 
iterative dynamics. The idea is that with the soft high T dynamics, the V; are allowed 
to short-cut through an interpolating (off-shell) space. As T slowly falls to zero, the V; 
are eventually forced on-shell, and finally a local minimum of the on-shell energy function 
E{si) crystallizes out. This way, one can hope to obtain better minima then by just using 
gradient descent. 



3 The Lennard-Jones Potential 

In models of proteins one often uses a potential between the individual atoms consisting 
in a sum of "bonded" interactions and "non-bonded" ones. The latter consist of Coulomb 
and LJ interactions between all the atoms. 

It has been argued [12] that a reasonable simplification results from considering the amino 
acids as effective monomers, interacting via effective potentials. We will focus on effective 
LJ potentials. Crucial physical properties like hydrophobicity then are efficiently incorpo- 
rated by a suitable choice of LJ parameters. 

The LJ potential between two monomers k and I is given by 



U; = /3v; - V;£(vi, . . . , Vjv) 



(13) 
(14) 




V kl (r) 




(15) 



An example is shown in fig. 3. It is short-range, and has a very strong (and unphysical) 
short-distance repulsion. If Aki > it has a local minimum located at 
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Figure 3: The Lennard-Jones potential (eq. (15)) and a regularized version (eq. (18)) as full and dashed 
lines, respectively (R=l; A=2). 

The energy to be minimized is given by the sum of all pair-potentials 



3.1 Regularizing the r — > Behaviour 

The very steep short-distance behaviour of the LJ-potential might cause problems for 
updating in minimization algorithms - not only for our variational approach but also for 
e.g. some Monte Carlo (MC) algorithms. A convenient way of controlling the short- 
distance singularity is by using a modified potential (see fig. 3) 



where is the minimum value of Vki{r). The modified potential only has a logarithmic 
singularity at the origin, and yields the same result in the limit 7 — > as V . Furthermore, 
the position of the minimum is preserved for all values of 7. The idea is to start with 
7 > and then gradually decrease 7 to 0. In minimization schemes of annealing type, 
where noise is introduced via a temperature T, the decrease in 7 should be coupled to the 
annealing in T (see Sect. 4). 

3.2 Coupling Distributions 

We choose the LJ parameters Rm and Am in two different ways; (A) identical values for 
all monomers, corresponding to dense packing, and (B) using empirical couplings for a 
random sequence of aminoacids, based on empirical frequencies of occurrence: 



(17) 




(18) 



7 



A For all pairs, set Rki = 1 and Aki = 2, leading to rj^ = 1. This choice implies an 
approximately dense-packed ground state. 

B Choose Rki and Aki in a phenomenological way as follows. Generate a random se- 
quence of amino-acids according to the known frequencies of occurrence in proteins 
[13]. The various values for Rki and Aki are derived from the empirical effective LJ 
couplings given in [12]. The resulting Aki distribution is shown in fig. 4. We note that 
while this distribution exhibits quite some width, it is far from being as dramatic as 
in the toy models proposed in refs. [8, 9], where Rki = 1 and Aki = 3.8 + 6.0 rand[0, 1] 
were used. 



0.00 




2.5 5 7.5 10 

r~ 6 coefficient x 10~ 6 
Figure 4: Distribution of Aki according to option B. 



4 The Variational Rotor Algorithm 



We next briefly describe how we implement eqs. (13,14) for the regularized LJ-potential, 
eqs. (15,18). The non-uniqueness due to eq. (11), and the existence of unique pairwise 
local minima rj^ in eq. (16), suggests a modified distance expression: 




vf 



?) (19) 
where a(kl) denotes the set of bonds connecting monomers k and I. To avoid instability 



ta- 



in the dynamics, rj^ is maximized to y \l — k\; this value is also used for a negative A 
On-shell, this yields the correct distance, while at the beginning of the simulation when 
V; = 0, most pairs are formally at their distance of minimum energy. The regulator 7 is 
chosen to depend on Rki as 

7ta = I 1 (20) 
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where 7 is the same for all pairs; this leads to correct relative normalization of the regu- 
larized potentials at short distances. 



The resulting variational equations (13,14) read 



u; = /3vi + 6 



\k-i\ 



Vi = g( Ui /T) 

where a denotes a connected subset of the bonds, that includes bond i. 
The variational algorithm then takes the following form: 

1. Initialize: 

1.1 Set T=T and 7=70- 

1.2 Initialize V; to small (~ .01) random vectors. 

2. Repeat until ^Ev- > 0.99999: 

2.1 Update all V;'s according to eqs. (21) and (22). 

2.2 Anneal: T = kxT;^ = fc 7 7 



(21) 
(22) 



The choice of parameter values are shown in table 1; for the stabilizer /3 the choice is 
motivated by the empirical observation that the /3-values required for high quality solutions 
scale linearly with N . In the gradient descent limit, where T = and the rotors are on-shell, 
V; = Uj/|uj|, j3 plays the role of an inverse step length. 





P 


7o 




To 


kx 


dense 


25N 


0.5 


0.99 


(3/3 


0.995 


real 


1007V 


0.5 


0.99 


/3/3 


0.995 



Table 1: Parameter settings for the variational algorithm (the first three apply also to GD). 



5 Numerical Results for the LJ-Potential 



We have compared the performance of the variational rotor algorithm to that of (i) a 
GD algorithm and (ii) a SA method - both for the model A, in which the final structure 
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will be close to a dense-packed configuration, and for the phenomenological model B. The 
GD method was chosen as the T = limit of the variational algorithm, with identical 
parameter settings for /3, 70 and fc 7 (cf. Table 1); thus, the regularized LJ potentials were 
used also here. The SA simulations were made using Metropolis updating, with an initial 
temperature T given by three times the modulus of the lowest energy achieved by the 
variational method in case A (for simplicity, the same T was used also for the B problems 
of the same size). With an annealing rate per sweep of kx = 0.99999, approximately 2.5 • 10 6 
sweeps were required, with a sweep defined by N attempted single-rotor updates. 

These choices of parameters and number of sweeps were based on reasonable trade-offs 
between solution quality and consumption of CPU time. 





CPU-time 


N = 10 


20 


40 


GD 


12 


70 


700 


Var. 


15 


120 


1500 


MC 


4000 


25000 


200000 



Table 2: Approximate CPU consumption for model B in seconds on a DEC Alpha workstation for the 
different approaches. Similar numbers hold for model A. 



The reference values for the ground-state energies thus obtained by SA, at the price of very 
high CPU consumption (see table 2), are only rarely equalled by the energies obtained using 
the two deterministic approaches. In model A (dense-packed) the variational algorithm 
clearly outperforms the GD one (see table 3). However, for the slightly more realistic 





N 


(Eqd — E var ) 


± 


(E var ) — EsA 


± 


ft probl. 


ft runs 


ft sweeps 
Var. GD 


A 


10 


0.8 


0.1 


1.0 


0.1 


1 


100 


2200 550 


20 


3.0 


0.4 


4.8 


0.2 


1 


100 


2300 1300 


40 


7.7 


0.8 


6.4 


0.3 


1 


100 


3000 4400 


B 


10 


0.6 


0.5 


6.0 


0.9 


10 


10 


4100 3800 


20 


2.3 


1.2 


14.9 


1.8 


10 


10 


9200 10400 


40 


9.1 


4.3 


59.7 


6.1 


5 


10 


16100 15300 



Table 3: Average performance differences. The number of runs refers to the variational and GD approaches. 
For each problem a single run was performed with SA. 



model B the two approaches seem comparable. A possible explanation for this is that 
the main advantage of the variational approach is the extra degrees of freedom since the 
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rotors can go off-shell; this facilitates the escape from local minima, which probably is 
more important in the dense packed case. 



6 Polyelectrolytes 



In this section we briefly discuss how the variational rotor approach can be applied to the 
case of a polyelectrolyte at a finite temperature. Theoretically, the approach suffers from 
a certain arbitrariness, but it gives quite good results for the case of unscreened Coulomb 
forces. 

The interaction energy for a polyelectrolyte consisting of N monomers with Coulombic 
repulsion forces is given by 

i? = £-i (23) 

kl T kl 

The variational average (E)v is then perfectly convergent. It is however difficult to eval- 
uate, and we will for that reason also here employ the simplifying trick of eq. (11). In 
this case, the resulting non-uniqueness is used to replace v? by 1 everywhere; thus, rki will 
evaluate to 

rki = J\k-l\ + J2 v i- v i ( 24 ) 

where i,j are restricted to the set of bonds linking monomers k and I. 

Then it is not a priori clear to what degree the results will make sense for a finite T . 
However, the replacement v? — > 1 gives the correct result for (t"|j)v, and we argue that 
it should be a fair approximation also for (l/rki)v — it is certainly correct for T — > 
(|vj| — > 1), and qualitatively correct for large N when T — > oo (v; — > 0). 

Thus armed with some confidence, we have used this approach to generate configurations 
in terms of {v;} for polyelectrolytes ranging in size from 20 to 1024, by iterating equations 
analogous to eqs. (21,22). We have chosen to characterize the configurations by their r.m.s. 
end-to-end distance r ee , which in the variational approach is given by 



rl 



(tf - 1) + 5>i ■ v,- (25) 



As can be seen from table 4, the results are in surprisingly good agreement with Monte 
Carlo (MC) data. 

We have also attempted a similar comparison for polymers with a screened Coulomb in- 
teraction. The results in that case are not nearly as good as in the unscreened case, and it 
is not difficult to see why. With screening, the interaction turns short-range, and for large 
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N 


Var. 


MC 


% dev 


20 


11.53 


12.06 


-4.6 


40 


26.80 


26.43 


1.3 


80 


59.05 


57.07 


3.5 


160 


125.8 


122.1 


3.0 


256 


207.6 


202.5 


2.5 


512 


429.5 


422.0 


1.4 


1024 


880.3 


870.2 


1.1 



Table 4: End-to-end distance r ee for polyelectrolytes. Comparison of variational results with MC data for 
different N for T = 0.837808 (room temperature). The last column lists relative deviations between the 
variational method and MC results in percentage. The errors in the MC data are approximately 0.1 %. 



molecules, essentially all the rotors will be identical in size and direction, V; = v, with v 
independent of N. Thus, we get for the end-to-end distance 

rl a (N _ d + (" -!)(*-») ,, (26) 

which goes like N 2 (unless v = 0, which happens for high enough T). Thus, r ee will scale 
like N, which is clearly unphysical for a short range potential. This failure is probably 
more due to the linear Ansatz for Ey being unsuitable for this problem, and not so much 
to the crudeness of the additional approximation of eq. (11). 



7 Summary 



We have developed a variational approach for finding approximate energy minima for 
proteins modelled by polymers with Lennard-Jones pair interactions between pointlike 
monomers. It has been numerically explored for two cases — dense-packed systems and 
more elongated ones. 

In the latter case, phenomenological pair potentials were used with random amino-acid 
sequences based on natural frequencies of occurrence. This leads to effective coupling 
distributions far more narrow than those commonly assumed in generic investigations of 
spectra and stability issues [9, 8]. 

For dense-packed systems the variational performance compares favourably with a gradient 
descent method with respect to solution quality, whereas for the more elongated ones the 
gain is very small. We interpret this difference as being due to the fact that the elongated 
chain provides less of an optimization challenge as compared to the dense packed ones. 
Hence there is a less pronounced difference between the methods. 
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The deterministic variational method fails to find in a consistent way the low energy states 
produced by a stochastic simulated annealing algorithm. However, the latter requires a 
factor 100 more in CPU consumption. 

As a by-product we have devised a regularization of the Lennard- Jones potential, where the 
short-distance steepness is gradually turned on in the updating process, thereby avoiding 
ill-behaved dynamics. Most deterministic (and some stochastic) methods will benefit from 
such a regularization. 

The variational method has been applied also to the case of polyelectrolytes at non-zero 
temperature, where the procedure is somewhat more ambiguous. For the case of unscreened 
repulsions the method yields good results as judged by MC, whereas it breaks down for 
the screened case. 
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